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Abstract 

<N 

The parareal in time algorithm allows to efficiently use parallel com- 
puting for the simulation of time-dependent problems. It is based on a 
\^ decomposition of the time interval into subintervals, and on a predictor- 

corrector strategy, where the propagations over each subinterval for the 
I corrector stage are concurrently performed on the processors. 

In this article, we are concerned with the long time integration of 
_^ Hamiltonian systems. Geometric, structure-preserving integrators are 

preferably employed for such systems because they show interesting nu- 
. merical properties, in particular excellent preservation of the total energy 

of the system. Using a symmetrization procedure and/or a (possibly also 
symmetric) projection step, we introduce here several variants of the orig- 
CS inal plain parareal in time algorithm [281 [2] [4] that are better adapted to 

the Hamiltonian context. These variants are compatible with the geomet- 
ric structure of the exact dynamics, and are easy to implement. 

Numerical tests on several model systems illustrate the remarkable 
properties of the proposed parareal integrators over long integration times. 
Some formal elements of understanding are also provided. 

Keywords: parallel integrators, Hamiltonian dynamics, long-time inte- 
gration, symmetric algorithms, symmetric projection, geometric integra- 
tion. 
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1 Introduction 



Increasingly intensive computations now become possible thanks to the improve- 
ment of both the efficiency and the clock rate of processors, the interprocessor's 
connections and the access to the different levels of memory. In addition, par- 
allel computing platforms, which allow many processors to work concurrently 
also become available. This second feature can only be useful if the problem to 
be solved can be decomposed into a series of independent tasks, each of them 
being assigned to one of the processors. 

The design of efficient algorithms for parallel architectures is the subject of 
intense current research. In the case of models governed by partial differential 
equations, most of — if not all — the contributions of the last three decades 
perform domain decomposition. We refer to [331 I40j for a review on recent 
advances, and also to the proceedings of the Domain Decomposition Method 
meetings (see www.ddm.org) for various achievements. When the problem is 
time-dependent, or when the problem is solely governed by a system of ordinary 
differential equations, relatively few contributions are available. We refer e.g. to 
the book of K. Burrage [9] for a synthetic approach on the subject (see also [TO]). 
In this book, the various techniques of parallel in time algorithms are classified 
into three categories: (i) parallelism across the system, (ii) parallelism across 
the method and (hi) parallelism across time. 

The parareal-in-time method, our focus here, was first proposed in the work 
of J.-L. Lions, Y. Maday and G. Turinici in 2001 28 . It belongs to the third 
category where parallelism is achieved by breaking up the integration interval 
into sub-intervals and concurrently solving over each sub-interval. The obvi- 
ous difficulty is to provide the correct initial value for the integration over each 
of these sub-intervals. Most of the techniques in the third category are mul- 
tishooting techniques, starting from the precursory work by A. Bcllen and M. 
Zennaro [TJ. This next led to the waveform relaxation methods introduced by 
E. Lelarasmee, A.E. Ruehli and A.L. Sangiovanni-Vincentelli [27], and to the 
multigrid approaches introduced by W. Hackbush [TS] . See also [TTJ . 

As it has been explained by M. Gander and S. Vandewalle in [TS], the parareal 
in time algorithm can be interpreted both as a type of multishooting technique 
and also as a type of multigrid approach, even though the bottom line of the 
approach is closer to that of spatial domain decomposition with coarse grid 
preconditioners . 

It is intuitively clear why so few contributions propose parallel-in-time algo- 
rithms: time-dependent problems are intrinsically sequential. On the other 
hand, the development of parallel computing provides computational oppor- 
tunities that exceed the needs of parallelization in space. This motivates the 
development of efficient parallel-in-time approaches. 

The parareal in time algorithm is based on the decomposition of the temporal 
domain into several temporal subdomains and the combination of a coarse and 
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a fine solution procedure. The name "parareal" has been adopted to indicate 
that the algorithm has, in principle, the potential to so efficiently speed up the 
simulation process that real time approximation of the solution of the problem 
becomes plausible. Since the original work |28j . where the convergence is proven 
to be linear, the parareal algorithm has received some attention and new de- 
velopments have been proposed. In ]3J, G. Bal and Y. Maday have provided 
a new interpretation of the scheme as a predictor-corrector algorithm (see also 
L. Baffico et al. [2]). The scheme involves a prediction step based on a coarse 
approximation for a global propagation of the system and a correction step com- 
puted in parallel and based on a fine approximation. G. Bal (in [3J) and E.M. 
Ronquist and G. Staff (in [35]) next provided some analysis on the convergence 
and the stability of the scheme. In [TS], M.J. Gander and S. Vandewalle proved a 
superlinear convergence of the algorithm when used on bounded time intervals, 
and a linear convergence on unbounded time intervals. 

In the past few years, the parareal algorithm has been successfully applied to 
various types of problems (see [55] for a review): nonlinear PDEs [3J, control 
problems [3D], quantum control problems [3T], Navier Stokes equations [T3] . 
structural dynamics [12j . reservoir simulation [17] . 

Although the plain parareal algorithm has proved to be efficient for many time- 
dependent problems, it also has some drawbacks in some specific cases. This 
is for instance the case for molecular dynamics simulations (as pointed out 
in [2]), or, more generally, for Hamiltonian problems. The exact flow of the 
system then enjoys many specific geometrical properties (symplecticity, possi- 
bly time-reversibility, . . . ). Some quantities are preserved along the trajectory: 
the Hamiltonian (e.g. the total energy of the system), and, in some cases, other 
quantities such as the angular momentum, . . . See the textbooks [2DJ [35] US] for 
a systematic introduction to numerical integration techniques for Hamiltonian 
systems. Symplectic and symmetric integrators are known to be suitable inte- 
grators for Hamiltonian systems. It turns out that, even in the case when the 
parareal algorithm is based on coarse and fine integrators that enjoy adequate 
geometrical properties (such as symplectic or symmetric integrators), the global 
parareal algorithm itself does not enjoy any of these properties. Consequently, 
the long time properties of the numerical flow (including e.g. energy preserva- 
tion) are not as good as expected. In this article, our aim is to design a parareal 
scheme that preserves some geometrical structure of the exact dynamics. 

Our article is organized as follows. In Section[2] we briefly recall the parareal al- 
gorithm, as presented in [2][3J. We next discuss, in Section[3j general, commonly 
used numerical schemes for Hamiltonian problems. We show the deficiencies of 
the parareal algorithm on such problems. In Section [4] we develop a symmet- 
ric version of the parareal algorithm, in a sense made precise at the end of 
Section |4.1[ As an alternative to symmetrization, we explore in Section [5j the 
idea of projecting the trajectory onto the constant energy manifold. Next, in 
Section [6] we couple a projection step with the symmetric algorithm developed 
in Section [4] while keeping the overall scheme symmetric. Combining these 
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two ideas, symmetry and projection, yields the most efficient algorithms. The 
performances of all the schemes proposed in these sections are illustrated by 
numerical simulations on two low-dimensional systems, namely the harmonic 
oscillator and the two-dimensional Kepler problem. The computational com- 
plexities of the different parallel integrators are analyzed in details in each of 
the respective sections. In Section [7J we consider a test case in a higher dimen- 
sional phase space, namely the simulation of the solar system (see also [35] for 
the derivation of algorithms specific to this test case, involving some parallel 
computations). We demonstrate that the efficiency and the qualitative prop- 
erties observed in the previous test cases carry over to this more challenging 
example. Using the symmetric parareal scheme with symmetric projection de- 
veloped in Section [6] (see Algorithm [3]) , we obtain a speed-up of more than 60 
with respect to a fully sequential computation (provided a sufficient number 
of processors are available) , for an equal accuracy. Section [8] summarizes our 
conclusions. 



2 The plain parareal algorithm 

In this section, we review the plain parareal method for a time-dependent prob- 
lem in a general setting. Consider u solution to the Cauchy problem 

du 

— + f(u) = 0, t € [0,T], supplied with the initial condition u(0) = u Q . (1) 

We assume standard appropriate conditions on / ensuring existence, uniqueness 
and stability with respect to perturbations, of the solution u to this problem. 
Let £ be the exact propagator, defined by £ t (uq) — u(r), where u(t) denotes 
the solution at time r of problem ([I]). 

In the sequel, for reference, we approximate the exact propagator £ by using an 
accurate numerical scheme. We can use any of the classical one-step schemes 
(explicit or implicit Euler schemes for simple problems, velocity Verlet scheme 
in the case of Hamiltonian dynamics) with a sufficiently small time step St. 
The associated discrete propagator is denoted as J 7 , with no reference to the 
time step St, in order not to make the notation too heavy. F t (uq) is thus 
an approximation of £ t (uq). Assuming that the approximation u n of u(T n ) is 
known at some time T„, the approximation of u(T n+ \) at a latter time T n +i 
is computed by performing (T n+ i — T n )/5t steps of the fine scheme, that is 
denoted, with the previous notations, 

U n +1 =^T 11 + i-T„(«n)- 

In the sequel, we will consider the dynamics 

q = M~ x p, p=-W(«), (q,p) eR d xR d , 
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where M is a diagonal matrix, and V is a smooth scalar function depending 
on q. We will integrate this dynamics using the velocity Verlet scheme, which 
is explicit and of order 2, and reads 

Qn+l 
Pn+1 

The plain parareal algorithm builds a sequence of iV-tuples u k = { u n} 1< „ <A r 
that converges, when k — ¥ oo, to the solution given by the fine scheme T: 
limfc-j.oo u\ = u n . In the sequel, we will consider Hamiltonian problems, for 
which designing schemes with a non-uniform discretization is not straightfor- 
ward. We thus restrict ourselves to the case of a regular discretization, namely 
T n = nAT with AT = T/N for some N £ N* , where [0,T] is the time interval 
of interest. We introduce another approximation Q of the exact flow E , which 
is not as accurate as J 7 , but is much less expensive to use than J- '. For example, 
we can choose the same discretization scheme as for J 7 , but with a larger time 
step dT 3> St (see Figure [T]). Hence, computing Q T amounts to performing r/dT 
time steps of length dT. Here again, the dependency of Q with respect to dT 
is not explicitly written. Another possibility is to define the solver Q from a 
simpler problem, which does not contain as much information as the original 
problem, and thus is easier to solve. We will use this opportunity in Section [7j 
In all the other sections, the only difference between Q and T lies in the choice 
of the time step. 

a 

i i i i 

dT 



In 



M~ 

St 



St p n 



se 



(2) 



Pn~ -(VV{q n )+VV(q n+ i)). 



Figure 1: Decomposition of the interval [0, T] in sub-intervals [T n ,r n +i], on 
which the solution can be integrated using a coarse scheme of time step dT, or 
using a fine scheme of time step St. 

Assume that we know an approximation {u^} Q<n<N of the solution to ([T]), at 
the end of some iteration k. Then the parareal scheme defines the next iteration 

{ U n +1 }o<n<N ^ 

u k n X\ = QA T {u k n +1 ) + Jat(^) - Sat(u£). (3) 

with the initial condition Uq +1 = uq. At the beginning of iteration fc + 1, we 
first compute .Fat(m^) — Qat(Uu) m parallel over each interval /„. Once this 
is completed, we only need to compute Gat{u^ v ) and add to it the stored 
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correction Fat(u^) — Gat{u*)- This is a sequential process, the complexity 
of which is negligible compared to the computation of Fat Note that 

an improved implementation in parallel has recently been proposed in [6] and 
explained hereafter. The analysis of the complexities of the different parareal 
schemes proposed in the sequel will use this implementation. 

From the construction of the scheme, convergence in at most N iterations can be 
demonstrated [29]. It can also be proven that, on a fixed time interval [0, T] and 
under some regularity conditions, the scheme ^ yields a numerical solution vfc 
after k iterations that approximates u(T n ) with an error which is upper bounded 
(up to a multiplicative constant independent of the time steps 8t, dT and AT) 
by err F + [err Q] k , where err F is the error on the solution at time T for the 
fine solver and err Q is the error on the solution at time T for the coarse solver 
(see [281 EES])- 



3 Parareal integration of Hamiltonian systems 

As announced in the introduction, the purpose of this article is to design parallel- 
in-time numerical schemes, derived from the parareal in time algorithm (|3|, for 
the integration of Hamiltonian dynamical systems. We first review the speci- 
ficities of such dynamics before discussing their numerical integration with the 
parareal algorithm 



3.1 Numerical integration of Hamiltonian systems 

In this article, we consider finite dimensional Hamiltonian systems, namely dy- 
namical systems that read 

• 9H . dH 

q= w p= ~w (4) 

where the Hamiltonian H(q,p) is a smooth scalar function depending on the 
variables q — (qi, . . . , qa) £ R d and p = (pi, ■ ■ ■ ,Pd) € 

The evolution of many physical systems can be written as a Hamiltonian dy- 
namics. Examples include systems in molecular dynamics (where q and p re- 
spectively represent the positions and momenta of the atoms composing the 
system, see [13]). celestial mechanics (where q and p represent the positions and 
momenta of planets or satellites [531 EH] ) , solid mechanics (after space dis- 
cretization, the wave equation modelling clastodynamics yields a Hamiltonian 
dynamics of the type Q, see [22112I]). -"- n a ^ these cases, H(q,p) is, physically, 
the total energy of the system. 

Hamiltonian dynamics have very specific properties that we now briefly review 
(see [20] [38] for more comprehensive expositions). First, the energy H(q,p) is 
preserved along the trajectory of Q. Second, the flow of a Hamiltonian system 
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is symplectic. It is well-known that symplectic schemes, such as the velocity 
Verlet scheme ([2]), are appropriate schemes for long time numerical integration 
of Hamiltonian dynamics. They indeed define a symplectic numerical flow, 
and thus preserve, at the discrete level, a fundamental geometrical property of 
the exact flow. In addition, the numerical flow given by a symplectic scheme 
applied on the dynamics Q can be shown to be almost equal to the exact 
flow of a Hamiltonian dynamics, for a so-called modified Hamiltonian, which is 
close to the original Hamiltonian H of Q . This is one of the main results of the 
celebrated backward error analysis for Hamiltonian dynamics (see (HUH [34], and 
the comprehensive survey in 20, Chap. IX]). As a consequence, preservation of 
the energy can be shown, in the sense that, for all n such that nSt < Cexp(c/<5t), 

\H(q n ,p n )-H(q ,p )\<C(StY, 

where 5t denotes the time step, r the order of the scheme, and the constants c 
and C are independent of St. The numerical error on the energy hence does not 
grow with the simulation interval (at least for time intervals exponentially long 
in the inverse of the time step). In some cases (namely, the case of completely 
integrable systems and of almost completely integrable systems, such as the 
solar system), the difference between the numerical and the exact trajectories 
can be shown to grow linearly in time, rather than exponentially, as would be 
the case for a generic integrator used on a generic dynamics u — f(u). 

In this article, we consider Hamiltonian systems for which the energy reads 

H(q, P ) = ^p T M- 1 p + V(q), (5) 

where, as above, M is a diagonal matrix, and V is a smooth scalar function 
depending on the positions q. For such an energy, the dynamics Q reads 

q = M- x p, p = -VV{q). (6) 

Setting u = (q,p), we recast this system as u = f(u). We then observe that 
the system is reversible with respect to p defined by p(q,p) = {q, —p), that is 
/ satisfies / o p = —p of. As a consequence, for any t, the exact flow ^ = £ t 
of Q is p-reversible, that is 

po$ = r 1 op. (7) 

Another class of interesting schemes for reversible Hamiltonian systems such 
as Q is the family of p-reversible schemes: a one-step scheme defined by 
u n+ \ — $>st( u n) is p-reversible if ^st satisfies Properties similar to the 
ones obtained with a symplectic scheme (good preservation of the energy, . . . ) 
can be shown when a reversible Hamiltonian dynamics is integrated using a 
p-reversible scheme. Conditions under which such properties are proven are 
however more restrictive than those of the backward error analysis for symplec- 
tic schemes: the Hamiltonian system should be reversible integrable (this implies 
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that, if q and p are in M. d , then the dynamics preserves at least d invariants) 
and a non resonant condition should be satisfied [2.0, Chap. XI]. 

The construction of p-reversible schemes is often done by first designing sym- 
metric schemes. A one-step scheme defined by u n +i = ^st( u n) is symmetric 
(see [101 Def. V.1.4]) if, for any St, 

t>st o = Id. (8) 

If a symmetric scheme satisfies the additional property 

)3°*St=*-i(Oft (9) 

then the scheme is p-reversible in the sense of ^ (see [SDJ Theo. V.1.5]). In 
practice, condition ^ is much less restrictive than the symmetry condition ([8]), 
and is satisfied by many schemes (including e.g. the forward Euler scheme). 
In the sequel, we will design symmetric schemes, satisfying and check a 
posteriori that they indeed satisfy Q. These schemes are hence p-reversible, 
namely satisfy Q. 

Despite the restrictions on the type of Hamiltonian systems for which p-reversible 
schemes allow for good long-time properties, symmetric schemes represent a very 
interesting alternative to symplectic schemes, because they are often easier to 
design. 

We conclude this section by recalling that the velocity Verlet scheme ^ , which 
is symplectic, as pointed out above, is also symmetric and p-reversible. 



3.2 Parareal integration of Hamiltonian dynamics 

We observe that, even if the coarse and fine propagators employed in the def- 
inition of the parareal algorithm ^ are symplectic or respectively symmetric, 
then, for a given parareal iteration k > 1 and at a given time step n, with 
n > k, the application uq = (qo,Po) l— ^ u n — (in^Pn) 1S neither symplectic nor 
respectively symmetric. 

This lack of structure has immediate practical consequences when employing 
the plain parareal algorithm Q on Hamiltonian systems. Consider as a first 
example the one-dimensional harmonic oscillator, with Hamiltonian 

H{q,p) = \p 2 + \q\ peR,qeR, (10) 

that we integrate up to time T = 10 4 . Set AT = 0.2, and consider the velocity 
Verlet scheme (which, we have recalled, is both symplectic and symmetric) for 
the fine and the coarse propagators, with time steps St — 10 -3 and dT = 0.1, 
respectively. 
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As a confirmation of the lack of geometric structure, we observe the lack of 
energy preservation. On Fig. [2] we plot the relative error on the energy preser- 
vation, defined, at time nAT and iteration k, as 

fc \H(qtPn)~ H(q ,Po)\ 

eTI n = J [777 \1 • ( n ) 

\H(q ,p )\ 

We note that, for short times (say up to T = 10 3 for the iteration k = 5), energy 
preservation is equally good for the parareal algorithm and for the sequential 
fine scheme. However, it deteriorates for larger times. Not unexpectedly, we 
do not observe the traditional behavior of geometric schemes (either symplcctic 
or symmetric) , namely a good preservation of energy, even when the numerical 
trajectory is very different from the exact trajectory. 

On Fig. [3j we plot, for several iterations k, the error 

erT k n = \\q k n -<ff\\ + \\p k n -p™ { \\ (12) 

on the trajectory (q,p), with respect to a reference trajectory (q re uPrei) com- 
puted with the velocity Verlet algorithm — used sequentially — with the time 
step St/10, where St is the time step of the fine propagator T . For k — 5, on the 
time interval [0, 10 3 ], we note that results are very good: the parareal trajectory 
is very close to the fine scheme trajectory, for a much smaller computational ef- 
fort. Indeed, at k = 5, assuming that the complexity of the coarse propagations 
is negligible, each processor has only computed 5 fine scale trajectories of length 
AT, in contrast to the situation when a sequential algorithm is used, and the 
processor has to compute T/ AT = 5000 fine scale trajectories of length AT to 
reach the time T = 10 3 . However, we also see that the error at iteration k = 5 
increases for t > 10 3 , and becomes unacceptable for times t of the order of 10 4 . 
With more iterations (say fc = 15), convergence of the trajectories up to the 
time T = 10 4 is obtained. 

Similar observations hold for the two-dimensional Kepler problem, where 

H(q, P ) = \p t p- i, P eR 2 ,qeR 2 . (13) 
We do not include them here for the sake of brevity. 



3.3 Evaluation of the complexity of the plain parareal al- 
gorithm 

In this section, for future reference, we assess the computational complexity of 
the parareal algorithm ([3| for the integration of the Hamiltonian dynamics 
We perform this evaluation under the assumptions that (i) the coarse and the 
fine propagators integrate the same dynamics, using an algorithm that requires 
the same fixed number of calls (set here to one) to the right hand side of ^ per 
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harmonic: dl=1 ,0e-3. dT=1 .De-2. DT=2.0e-1 , PS 




harmonic: dt=1 .Oe-4, dT=1 .De-2, DT=2.0e-1 , PS 



k=1 1 
k=13 
k=15 




Figure 2: Error ( 11 ) on the energy for the harmonic oscillator problem, obtained 
by the parareal method ^ with St = 10" 3 , dT = 0.1, AT = 0.2. 
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harmonic: dt=1 .Oe-3, dT=1 .Oe-2, DT=2.0e-1 , PS 




time step (we consider that the complexities of evaluating \7H or W are equal, 
and we denote it by Cy), and that (ii) we have all the necessary processors to 
perform all the fine scale computations in parallel. 

We use here the implementation quoted to us by [B]. It consists in starting the 
computation of Jat(«' +1 ) immediately after u^ +1 has been computed in 
and not waiting that all (M* +1 )i< n <jv are available. This allows to start the 
parareal iteration k + 2 much sooner. The complexity of the first coarse propa- 
gation scales as Cy Tj dT. We next distinguish two extreme cases, whether the 
complexity of the fine propagator on a time interval of size AT is smaller (re- 
spectively larger) than the complexity of the coarse solver on the whole interval 

AT T 

[0, Tl. The first case corresponds to the situation when — - — < — , the second 

St dT 

AT T 
case when — — > — . 

St dT 

In the first case, the complexity of each iteration is dominated by the complexity 
of the coarse solver on the whole interval [0,T]. The complexity of such a 
parareal iteration is thus again of the order of Cy T/dT. In the second case, 
the complexity of each parareal iteration is of the order of Cy AT /St. In both 
cases, a final coarse iteration needs to be performed. 

Denoting by Kp the number of parareal iterations, the approximate complexity 
of the scheme ^ is of the order of 

Cp = (if P + l)Cy 1 



in the first case, and 



in the second one. 



dT 
AT 



Cp = Kp Cy 



In comparison, the complexity of the fully sequential algorithm, using the small 
time step St, is 

T 

C scq = Cy -. 

Note that the second case above corresponds more to the paradigm of the 
parareal scheme (remind also that, if possible, the coarse solver should be based 
on a less expensive dynamics than the fine solver, as in Section [7] below), espe- 
cially for the integration of large Hamiltonian systems as the one considered at 

T 

the end of the article. In that case, the speed-up is of the order — — , which 

KpAT 

is the number of processors divided by the number of iterations. 

In the sequel, for the complexity analysis, we shall assume that we are in the 
second case above, i.e. we assume that 

~sT > dr- (14) 
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4 A symmetric variant of the parareal in time 
algorithm 



In [5] , a symplectic variant of the parareal algorithm has been introduced. The 
strategy there is based on the reconstruction of the generating function S : 
(q,p) £ M d x R d h> K associated with the symplectic map Fat ° Sat- ^n 
interpolation procedure is used to obtain an approximation of that generating 
function. The optimal way to perform the interpolation is still an open question, 
especially when working with high dimensional problems (d> 1). In this article, 
we limit ourselves to the design of a symmetric scheme, using the fact that it 
is often simpler to symmetrize a given scheme than to make it symplectic. We 
hope to return to symplectic variants in future publications. 



4.1 Derivation of the scheme 



Our idea is based on the following well-known observation. Consider a general 
one-step scheme £4+i = ^ AT{U n )- Then the scheme 

£4+1 = *AT/2 ° (f-AT^r 1 (£/») (15) 

is symmetric (see [201 Chap. V]). For future use, we introduce the intermediate 
variables £7„+i/2 = f 1 ^— AT/2) 1 (£4), and write the above scheme as 

U n = *-AT/ 2 (^n+l/ 2 ), £4+1 = $AT/ 2 (^+l/2)- (16) 



Let us now write the parareal algorithm ([3| in a form more appropriate for 
our specific purpose. We first choose a number K of iterations for the parareal 
algorithm, and set 

Un ■= {Un,U l n , ■ ■ ■ ,U%). 

Then the parareal scheme ^ can be written as 

U n+ i = #at(C4), 



where the map Wat is defined by 

*Ar(Cn) - 



/ GatK) 

GAT(u n ) + Tat«) - GatK) 



V GAT(<)+T A T«- 1 )-gAT(u^ 1 ) j 



We now apply the symmetrization procedure ( |16[ ) to the map ^at, and consider 
the scheme 

U n = W-AT^ (£4+1/2), 



U n . 



AT/2 



(U n 



-1/2J 



(17) 
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Let us write this scheme in a more detailed manner. The first line of (17) reads 
..0 



Q-AT/2(ui +1/2 ), 
n = S-AT/2{u n+1 / 2 ) + -7 r -AT/2(M° + i/ 2 ) ~ G-AT/2 0" + 1/2): 



,K-1 



- Q- AT/2 (Un+ 1/2) + F-AT/2(u* + y 2 ) - G-AT/l{ 



K-2 



L K - 2 ) 
Vfl/2/' 

Un = G-AT/2(u* +1/2 ) + J r -AT/2{Un+l/ 2 ) ~ G-AT/2 (^+l/ 2 ) ■ 

We hence obtain 



-1/2 



*n+l/2 



l n+l/2 
i n+l/2 



S-AT/ 2 ( U ™)' 



-AT/2' 
'-AT/2 



— G-AT/2 



4 - J 7 -AT/2(u a n+1/2 ) + ^-AT/2(w° +1/2 ) 



- -? r -AT/2(«f +1 2 /2 ) + £-AT/2« +1 %) 



K-2 



'-AT/2 



,K 



-7 r -AT/ 2 (^ +1 1 /2 ) + e?-AT/2« + i/ 2 ) 



We now collect the above set of relations with the second line of ( 17 ) and obtain 
the following formulation: 



M n+l/2 — £-AT/2( M ")' U n+1 ~ QaT/2 ( U n +1/2 ) > 

at iteration k — 0, and, for k > 1, 



(18) 



Vfl/2 
"n+1 



— ^-AT/2 



- -7 r -AT/2« +1/2 ) + 0-AT/2« +1/2 ) 

Sat/2 «£/a ) +^at/2« +1/2 ) - £at/ 2 « +1/2 )- 



(19) 



In the sequel, the scheme (|18[)-( 19 ) is called the symmetric parareal scheme. 



Remark 1 We have recalled in Section 3.1 that p-reversible schemes have good 
geometrical properties. By construction, the scheme (l£fy -( 19) is symmetric. We 
have checked that it also satisfies condition ([9|. As a consequence (see again 
Section 3.1), the scheme | 18)- (Tify is p-reversible. 



We note that, apart from the computation of G_& T / 2 , the scheme ( 18 1-( 19 ) is 
completely explicit, as soon as the propagators Q and J- themselves are explicit. 
In particular, we point out that the inverse of the fine propagator J- is not 
needed. 

Let us now show that, in this new algorithm, all the expensive computations 
(involving J 7 ) can actually be performed in parallel. Assume that, at a given 

iteration k, we know |wj;i n ^ ^ AT and \ , , /0 > . Then the correction 

n) 0<n<N I ™+ 1 / 2 Jo<„<iV-l 

terms 

■?AT/2(\ + l/ 2 ) - QAT/2{u k n+ i/ 2 ) 
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and 

J 7 - AT/2 Kj+1/2 ) - 6-AT/2(Vl/2) 

can be computed independently, over each processor. The algorithm can next 
proceed with sequenti al, b ut inexpensive, computations. Note again that, as 
pointed out in Section 



3.3 



we can start the computation of J r AT/2(^ 



-1/2 



) and 



F-at/2(u, 1+1/2 ) as soon as u 



+ i/2 i s available (we do not have to wait for all 
the computations of iteration k to be completed). 



The symmetric parareal algorithm is summarized as follows. 
Algorithm 1 Let uq be the initial condition. 

1. Initialization: set Uq — Uq, and sequentially compute |w^ +1 y 2 | 



<n<JV-l 



and {<+i} 0<n<N _ 1 as 



-1/2 



"n+1 



GAT/2{u n+1/2 )- 



and{u* +1/2 } 

{ U nX\/2} 



0<n<N-l 
<n<N-l 



2. Assume that, for some k > 0, the sequences { M n}o< n <Ar 

are known. Compute the sequences { M n +1 } < n <jv an< ^ 
at iteration k + 1 by the following steps: 

(a) For all < n < N—l, compute in parallelJ-^AT/2( u n+i/2)' ^"AT/2( u „_|-i/ 2 ) 
£-AT/2«+i /2 ) and Qat/2{uI +1/2 ); 

(b) Set u* +1 = u ; 

(c) Compute, for < n < N — 1, 



,fe+i 

ri+1/2 
u n+l 



~ G-AT/2 



K + ~ ^-AT/2« + i /2 ) + 0-AT/2«+i /2 ) 
0AT/2(« n +l/ 2 ) +^AT/2K fe l+1/2 ) - QAT/2{u k n+1/2 )- 



Some comments are in order. The parareal scheme ^ is no t a o ne-s tep scheme 
defining v%\\ from u^ +1 , and hence the symmetric form (18l-(19l cannot be 
considered, strictly speaking, as a symmetric integrator in the classical sense. 
However, several reasons lead us to believe that this algorithm is the appro- 
priate generalization of symmetric integrators when dealing with parareal-type 
algorithms. 



First, the scheme at iteration k = 0, defined by (18 1, is exactly the symmetriza- 
tion (15 1 of the coarse propagator Gat- 



Second, the flow ( 18 )-( 19 ) is symmetric in the sense that, if for some n and 
some k, 

( U n+1> ' " ) u n+li u n+l) 



15 



is obtained from 



(u Q ■■■ u k u k+1 ) 



by the flow implicitly defined in ( 18 )-( 19 1, then 



(u° ••• u k u k+1 ) 



is obtained from 



i u n+li ' 



*n+l) a n+l) 



by the exact same algorithm reversing the time. In fact, only the considera- 
tion (and the storage) of the last two iterates (in terms of "parareal iterates") 
(u k l ,u^ l +1 ) is required to perform the iterations, and the above argument. 

Third, if the coarse propagator happens to be identical to the fine one (which is 
of course not supposed to be the case for efficiency of the parareal integrator!), 



then the symmetrized form ( 19 ) reads 

fe+l _ n-l ( y k+l\ 
"n+l/2 ~~ y_AT/2V u n )l 



,fe+l 
l n+l 



SaT/2{' 



,k+l 
n+l/2 



), 



and it thus coincides with the standard symmetrized version of the coarse prop- 
agator. 



Finally, formally taking the limit k 



-Si+l 



AT j2 



oo in (19| yields 

AT/2)- 1 ^). 



This shows that the limit of the symmetric parareal algorithm in terms of 
parareal iterations coincides with a standard symmetrized form of the map 
J- at- Note also that this latter algorithm is not the symmetrized version of the 
fine propagator, since symmetrization does not occur after each time step, but 
after AT/(2St) time steps. 

These observations show the formal consistency of our notion of symmetrization 
for parareal-type integrators with the classical notion of symmetry. 



Remark 2 Instead of considering (15) to symmetrize a given scheme "J 7 at, one 



can alternatively consider the symmetric scheme 

U n+1 = (*_at/ 2 ) _1 ° *AT/ 2 (!7n). 

In the parareal context, this yields an algorithm with similar properties as the 
above algorithm { 18)-\l(fy, and with comparable numerical results. In the sequel, 



for the sake of brevity, we only consider the symmetrization procedure (151 



4.2 Evaluation of the complexity of Algorithm [T] 

We assess the computational complexity of Algorithm[l]under the same assump- 
tions as in Section 3.3 We again assume (14), and we again start to compute 
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•^at/2(w*+i/2) anc * - 7r -AT/2( u ^ + i/ 2 ) as s00n as possible, e.g. as soon as u^ i+1/2 
has been computed. 

The complexity of the first coarse propagation scales as CV T/dT. The com- 
plexity of each parareal iteration is of the order of CV AT /St (note that we 
need to call twice the fine solver, but only for time intervals of length AT/2). 
A final coarse iteration needs to be done, the complexity of which is neglected 
in comparison to the cost of a fine propagation. Denoting by K$p the number 
of parareal iterations, we obtain that the complexity of this scheme is 

AT 

Csp = K SP Cy ~^r- 

The complexity is the same as that of the plain parareal algorithm. The speed- 
T 

up, which is equal to — — — — , is again equal to the number of processors divided 
a sp AT 

by the number of iterations. 



4.3 Numerical examples 



Let us now again consider the example of the harmonic oscillator ( 10 1, and into 



grate this Hamiltonian system with our newly constructed symmetric parareal 



algorithm (18)-(19). Note that system (10) is an integrable reversible Hamil- 
tonian system, and thus belongs to a class of problems for which symmetric 
integrators have been shown to preserve energy [20l Chap. XI], as has been 
recalled in Section [3711 



The relative errors (11 1 on the energy preservation are shown on Fig. |4j while 
the errors (12 1 on the trajectories are shown on Fig. [5j Parameters have been 
chosen such that the complexity to reach a given time step n at a given iter- 
ation k is equal for the symmetric parareal algorithm considered here and the 
standard parareal algorithm ^ discussed in the previous sections. Results are 
disappointing: they are indeed very similar, both qualitatively and quantita- 
tively, to those obtained with the plain parareal scheme ^ (see Section 3.2 
Figs. |] and [§. 

Similar results have been obtained for the Kepler problem, which is also an 
integrable reversible Hamiltonian system. 

It is useful and instructive to now explain such a poor behavior, using the 
following formal elements of analysis. We first observe that a parareal integrator 
may be seen, at parareal iteration k, as an integrator of a system consisting of 



k + 1 identical replicas of the original system under consideration. From ( 18 ) , we 
know that the first replica is integrated by a symmetric algorithm. If the system 
under study is an integrable reversible Hamiltonian system, its energy is thus 
preserved in the long time by the simulation at parareal iteration k — (this is 
confirmed by numerical experiments, see e.g. Fig. |4j curve k = 0). Next, since 
the replicas are noninteracting, the system of replicas is evidently an integrable 
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harmonic: dt=1 .0.6-3, dT=1 ,0e-2, DT=2.0e-1 , SPS 




harmonic: dl=1 ,0e-4. dT=1 ,0e-2. DT=2.0e-1 , SPS 




Figure 4: Errors on the energy for the harmonic oscillator problem, obtained 
by Algorithm [l] with St = 1CT 3 , dT = 0.1, AT = 0.2. 
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reversible system, and has an energy that is equal to k + 1 times the energy 



of the original system. Assume now that the symmetric propagator ( 18 )- ( 19 1 



applied to the system of k + 1 identical replicas conserves the energy of the global 
system (i.e. the sum Y!,e=Q H(q e ,p e ) of the energies of the replicated systems), 
as we could expect from a symmetric scheme applied to an integrable reversible 
Hamiltonian system. Under this assumption, it follows, by induction on k 7 that 
the energy is preserved along the trajectory, at each parareal iteration k. This 
is clearly in contradiction with the observed numerical results! The flaw in the 
above argument is that, for a symmetric scheme to preserve the energy of an 
integrable reversible system, we have recalled that, among other conditions, the 
frequencies present in the system have to satisfy a non-resonant diophantinc 
condition. This is precisely not the case here for the system of replicas, by 
replication of the original frequencies! 

The theoretical argument showing energy preservation does not apply, and nu- 
merical results confirm that energy is indeed not well-preserved, in the long-time 



limit. The symmetric parareal scheme ( 18 )- ( 19 1 hence needs to be somehow 
amended. 



4.4 Symmetric parareal algorithm with frequency pertur- 
bation 

To prevent the different replicas from being resonant, a possibility is to consider, 
at each parareal iteration fc, a system slightly different from the original system. 
For instance, in the case of the harmonic oscillator, we may consider at each 
iteration k a harmonic oscillator with a specific frequency u>^: 

Hk(q,p) = ^P 2 + ^lq 2 - 

The unperturbed case corresponds to ujk — inexact = 1 in the above energy. 
Provided the cok are all different from one another (and non-resonant), the 
system of replicas is non-resonant. In the following, the shift is chosen such 
that it vanishes when k — ¥ oo, i.e. lim^oo u>k = w eX act- 

More precisely, we introduce a fine propagator and a coarse propagator 

(k) 

f° r ^ ne Hamiltonian dynamics associated to the Hamiltonian Hk- Rather 
than symmetrizing the parareal algorithm ([3]), we consider the scheme 

77° - 

a n+l — »AT\"n/i 

_ G {k+1) (n k + 1 ) + T {k+1) (v k ') -G {k+1) (i, k \ [ ' 

U n + l — »AT \ U n I < ■'AT \ U nl » AT \ U n)-> 

(note the upper indices (fc + 1)) and symmetrize this scheme along the lines of 
Section |4~T| We thus obtain, at iteration k = 0, the algorithm 

"n+l/2 = (^AT/2) "n+1 -^AT/2( U n+l/2)' (^1) 
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and, for k > 1, 



,fc+l 



II 



1/2 

fc+1 
n+1 



= g 



,(fc+i) 

-AT/2 



,k + l 



-r 



(fc+i) 

-AT/2 



(«*, 



1/2-' 



r>(k+i) f 7 ,fc ^ 

y -AT/2^n+l/2' 



C ( ' £+1) (7/ fc + 1 1 , Tr( fc + 1 )|' 1 .fe \ _ C (fe+1) ('7/' £ 1 

y AT/2 \ a n+l/2> ~ ■ r AT/2 ^"n+1/2^ y AT/2 \ a n+l/2> 



(22) 



AT / 



'AT/2 ' 



Let us briefly discuss the consistency of this scheme. First, at iteration k + 1, 
if the coarse propagator is identical to the fine propagator F<- k+1 \ then 



(22) reads 



,k+l 
n+l/2 



-(k+1) 
-AT/2 



« 



k + l\ 



-1 _ tA»+V( u 
-1 — • r AT/2 \ u 



AT/2 \ u, n+l/2h 



which is a scheme consistent with the Hamiltonian dynamics driven by H^+i- 
Since the perturbation can be neglected when k — > oo, we recover, in the limit 
of large k, a scheme consistent with the original dynamics. 



In addition, formally taking the limit k — > +oo in ( 22 ) yields 



^AT/2°(^AT/2) 'K). 



where J 7 is the fine propagator associated to the original dynamics. This shows 



that the limit of the algorithm (21)- (22) in terms of parareal iterations coin 



cides with a standard symmetrized form of the fine propagator of the original 
dynamics. 

On Fig. [6j wc plot the errors (11) on the energy and the errors (12) on the 
trajectory, obtained with the algorithm (21 )-(22) applied to the harmonic oscil- 



lator problem. We focus on the results obtained at the last parareal iteration, 
the only one for which the perturbation vanishes. We observe that the en- 
ergy does not drift, while the trajectory is not extremely accurate (actually, 
both observations hold also for intermediate values of k). We are hence in the 
regime of geometric integration, with a good accuracy on energy being not a 
consequence of a good accuracy on the trajectory. This is an improvement in 
comparison to the previous algorithms (parareal algorithm ^ and symmetric 
parareal algorithm ( 18 )-( 19 1). 



We next turn to the Kepler problem, where we introduce a perturbation by 
considering 



U I \ 1 T a k 
Hk(q, P ) = 2 P 



(23) 



where the unperturbed case corresponds to = a oxac t = 1. Results are shown 
on Fig. [7] The same conclusions as for the harmonic oscillator test case hold. 

We hence numerically observe, in both cases considered, that the energy does 



not drift when we use the scheme (21)-(22), and that energy preservation does 



not owe to trajectory accuracy. This formally validates our understanding: the 



energy drifts observed with the symmetric parareal algorithm (18)-(19) are due 
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harmonic: dt=1.0e-3, dT=1.0e-1, DT=2.0e-1 , SPS + Shift(1.1, 0.9, 1.05, 0.95, 1) F shifting 



k=0 — 


k=2 ----- 








k=3 
k=4 




coarse 
fine 




i!fil'i™'-| 



















harmonic: dl=1 .Oe-3, dT=1 .0e-1 , DT=2.0e-1 , SPS + Shift(1 .1 , 0.9. 1.05. 0.95, 1) F shifting 




Figure 6: Errors on the energy (left) and errors on the trajectory (right) for the 
'lem, obtainei 
'turbation. 1 
0.95 and w fc=4 = w eX act 



harmonic oscillator problem, obtained with the symmetric parareal method (|21 1 

re u> k . 
1 [St 



(22J), with frequency perturbation. The frequencies are u>k=o — 1.1, Wfc=i = 0.9, 

10- 3 , dT = 0.1, 



^fc=2 = 1-05, w fe=3 
AT = 0.2). 
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Figure 7: Error on the energy (left) and errors on the trajectory (right) for 
Kepler problem, obtained with the symmetric parareal method (21 1-([22|, with 
frequency perturbation. The coefficient a k in (23 1 is cik=o = 0.9, otu=i = 0.95, 
dk=2 = 0.975, a k =3 = 0.9875, a k=i = 0.99375 and a fc=5 = a OX act = 1 (8t = 
IQ- 4 , dT = 10" 2 , AT = 0.2). 
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to resonances. As soon as the replicas are made non-resonant, we recover a 
behavior in terms of energy preservation that is typical of geometric algorithms. 

Nevertheless, even though energy does not drift, we also see that it is not ex- 
tremely well preserved. It is often the case that, at the final iteration, energy 
preservation is actually not really better than if obtained using the coarse propa- 
gator on the original dynamics. We also observe that the error on the trajectory 
remains quite large, for the practical values of k that we consider. So the algo- 
rithm ( 21 )-( 22 ) , although better than ^ and (18|-(19|, is not satisfactory. 



Remark 3 Another possibility to integrate at each parareal iteration systems 
that are non resonant from one another is to use, at each iteration k, a different 
time step dT^ for the coarse propagator. This strategy, that is more general 
than the one discussed above (indeed, "shifting the frequency" may not be easy 
if the system is not explicitly an harmonic oscillator), yields results that are 
qualitatively similar to the ones reported here. 

Remark 4 The idea of "shifting the frequency" is independent from the fact 
that the algorithm is symmetric. It is hence possible to apply the standard 
parareal method ^ on systems that are slightly different from one another at 



each parareal iteration (in (22), we have used the symmetric version of the 
parareal algorithm on such shifted systems). Results are slightly better when this 
strategy is used with the symmetric version of the algorithm. 



5 Parareal algorithm with projection 

As an alternative to symmetrization described in Section [4j we consider in this 
section a different idea. We couple the plain parareal method with a projection 
of the trajectory on a specific manifold, defined by the preservation of some 
invariants of the system (namely the energy, and possibly other ones). 

Observe that many systems have actually more than one preserved quantity: 
this is the case for the two-dimensional Kepler problem, where the angular 
momentum L = q x p and the Runge-Lenz-Pauli vector A = p x L — -r^-rr are 

h\\ 

also preserved, or for any completely integrable system. However, the generic 
situation is that some invariant quantities are known (among which the energy) , 
and that there may be other quantities that are preserved by the dynamics, but 
may not have been identified. As a consequence, we first consider a generic 
method, based on projecting on the constant energy manifold 

M = {(g,p)eR 2d ; H(q,p) = H a }, 

for the Kepler problem. Note that the harmonic oscillator is not a discriminating 
test case in this respect, since the energy is then the unique invariant. Next, for 
the sake of completeness, again for the Kepler problem, we consider a method 
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based on projecting on the manifold where both H and L are kept constant (see 
Section pTI ). 



Remark 5 Symplectic algorithms, such as SHAKE JgEjj or RATTLE fj Hffi . 
have been proposed for the numerical integration of Hamiltonian systems re- 
stricted to some manifold defined by holonomic constraints. Note however that 
here, the constraint at hand, H(q,p) = Hq, is not holonomic. 

Following [201 Chap. IV], we consider the following definition of the projection 
■km onto the manifold M. For any y = (q,p) £ R 2d , we define (q,p) = y = 
n M (y)£R 2d by 

V = n M (y) =V + XVH(y), where A G K is such that H(y) = H . (24) 
5.1 Parareal algorithm with projection 

The coarse propagator being the velocity Verlet algorithm, which yields an 
acceptable energy preservation, we do not need to project on the constant energy 
manifold at iteration k = 0. But, as shown by our numerical results of the 
previous sections, the trajectory departs from this manifold for the iterations 
fc > 1, in the long time limit. This motivates the consideration of the following 
algorithm. 

At iteration k — 0, we set 

«£?i=0at(u£=°). (25) 

For the next iterations, we simply couple the parareal algorithm ([3| with a 
projection step, as schematically represented on Fig. [8] We hence define the 
solution at iteration k + 1 from the solution at iteration k by 

<+! = km [QAT(u k n +1 ) + 7at(«J) - 0at(u*)] . (26) 




Figure 8: Schematic representation of the standard projection. 
This leads to the following algorithm. 

25 



Algorithm 2 Let uq be the initial condition. 



1. Initialization: set Uq = Uq, and sequentially compute { u n}i< <jv 
<+i = 5at«). 

2. Assume that, for some k > 0, the sequence { u n} 0<n<N * s known. Com- 
pute { u n +1 }g<„<jv °^ iteration k + 1 by the following steps: 

(a) For all < n < N — 1, compute Jat("b) — Gat(u^) in parallel; 

(b) Set u k +1 = u ; 

(c) Compute {un\\} a<n<N -i sequentially by 

<+! = *M [5at« +1 ) +Jat(<) - WO] ■ 

The additional complexity of the above algorithm, in comparison to the parareal 
algorithm ([3]), comes from solving the nonlinear projection step (24). Note how- 
ever that we expect to have a good initial guess for this problem, since the solu- 
tion at the previous parareal iteration is expected to be a coarse approximation 
of the exact solution. We hence expect to solve the nonlinear projection step 
with a few iterations of, say, a Newton algorithm. We will see in the sequel that 
this is indeed the case. 



5.2 Evaluation of the complexity of Algorithm [2] 

We assess here the complexity of Algorithm [2] under the same assumptions as 
in Sections |3.3| and |4.2| Again, the complexity of the first coarse propagation 
scales as C v T/dT. 

Let us denote by C pro j the complexity of each projection step. It depends on 
the number m pro j of iterations needed to solve this nonlinear problem. Each 
iteration requires to evaluate the gradient of the energy, a task the complexity 
of which is Cv- Hence, 

Cproj — Cv ^proj ■ 

In all the simulations we have performed, m pro j is rather small. This extra cost 
Cp ro j can be neglected with respect to the cost of the fine propagation. 

The complexity of each parareal iteration is thus again of the order of Cv AT /St. 
Denoting by K P / pm ^ the number of parareal iterations, we obtain that the com- 
plexity of this scheme is of the order of 

AT 

Cp/p X oj = K P/pIoj C V -jjr, 

as for the parareal scheme (l3l and its symmetric variant ( 18 )-( 19 1. The speed-up 

T . 
— — is again equal to the number of processors divided by the number 

A'p/projAT 

of iterations, as long as m pro j remains small. 
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5.3 Numerical results 



We consider the Kepler problem associated to the energy (13), and we plot 



the relative errors (11) on the energy preservation and the errors (12) on the 
trajectory on Figs. |9| and [T0| respectively 



The projection step (24) is solved using a Newton algorithm. Note yet that 
the energy is not exactly preserved, because only a few steps of the Newton 
algorithm are performed. Indeed, the algorithm is terminated as soon as one of 
the following three convergence criteria is satisfied: 



CI: the quantity err is smaller than the tolerance tol, 

C2: the maximum number N™£* of iterations has been reached, 

C3: the error on the energy does not decrease, 



where err is here the relative error on the energy ( 11 ). We have chosen to work 
with iV^ x = 2. It turns out that the best choice for tol is the relative error on 
the energy preservation of the fine scheme. In the current case, with St = 10~ 4 , 
this corresponds to taking tol = 10~ 7 . 

We observe an excellent energy preservation, although we possibly stop the 
Newton algorithm before convergence: the prescribed tolerance is reached over 
the complete time range for all parareal iterations k > 7. This is confirmed 
m Table [TJ where we show that, in the Newton algorithm, the criteria CI (the 
requested tolerance has been reached) is satisfied before the criterias C2 and C3. 



In addition, for k > 7, the error (27) on the angular momentum preservation 
is always smaller than 10 -2 , and it decreases down to 10~ 4 for the iterations 
k > 11. 



Criteria CI 


Criteria C2 


Criteria C3 


91.6 % 


6.8 % 


1.6 % 



Table 1: To stop the projection procedure, we use three criteria. We gather here 
how often each criteria triggered the projection procedure to stop in Algorithm[2j 



We also observe that only 11 iterations are needed to obtain convergence of the 
trajectory (that is, the trajectory obtained at iteration k = 11 is as accurate as 
the one given by the fine scheme) on the time range [0, 10 4 ]. 

With this algorithm, we are hence able to recover accurately the trajectory on 
the complete time range [0, 10 4 ], for a moderate number of iterations. Results 
are much better than with the other modifications of the parareal algorithm 
that we have described so far (for the same value of the parameters St, dT and 
AT). 
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*epier dt l Oe-4. dT i Oe-2. DT-2 Oe- 1 . PS +■ Projn 
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;ure 9: Errors on the energy for Kepler problem, obtained by Algorithm 
= 10" 4 , dT = 0.01, AT = 0.2). 
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On Fig. 11 we plot the relative errors on the angular momentum L(q, p) = qxp, 
which is another invariant of the Kepler problem. The relative error is defined 
by 

k _ \L(q k n ,p k n )-L(q 07 p )\ 

m%,po)\ 

This quantity does not blow-up. We yet observe that it is not preserved as well 
as would be the case with a geometric integrator. The rather good behavior of 
this invariant is a consequence of the accuracy of the trajectory. 



kepler: dl=1 .Oe-4, dT=1 .Oe-2, DT=2.09-1 , PS + ProjH 





Figure 11: Errors on the angular momentum for Kepler problem, obtained by 
Algorithm [2] (St = lfT 4 , dT = 0.01, AT = 0.2). 



Remark 6 Note that the parameters chosen for the numerical simulation of the 
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harmonic oscillator and the Kepler system do not satisfy assumption | H). Our 
aim with these two toy-problems is indeed more to understand the basic features 
of the variants of the parareal algorithm than to be as efficient as possible. Our 
viewpoint will be different when we consider a system of higher dimensionality, 
in Section^ below. 



5.4 Considering more than one invariant 

In this section, we again consider the case of the Kepler problem, and now take 
into account that, for this specific system, we know another invariant besides 
the energy, namely the angular momentum L = q x p. Define the manifold 

M 2 = {(q,p) € R 2d : H(q,p) = H , L(q,p) = L Q } , 

and consider the projection ttm 2 onto that manifold, defined by 

y = KM 2 (y) = v + AiVff (50 + A 2 vi(y), 

where (Aj., A2) <E K 2 is such that H(y) — Ho and L(y) = Lq. 

We couple this projection step with the parareal algorithm as follows . We first 



define tc at iteration k = by ( 25 ) 



u n+1 = y&T(u n ). 

We next define the solution at iteration k + 1 from the solution at iteration k 

by 

<X\ = *M, [Sat(^ +1 ) + J=AT«) ~ SMO] ■ (28) 
The algorithm is similar to Algorithm [2] where ttm is replaced by ttm 2 - 

Let us now integrate the Kepler dynamics with this algorithm (we again use 
a Newton algorithm to compute the Lagrange multipliers, with the same three 



stopping criteria as in the previous section). We plot the relative errors (111 



on the energy preservation, the relative errors (27) on the angular momentum 

re- 



preservation and the errors ( 12 ) on the trajectory, on Figs. 12 13 and 14 
spectively. 

We observe qualitatively the same behavior as for Algorithm [2] which is based 
on projection on the constant energy manifold. Only 8 iterations are needed 
to obtain convergence of the trajectory, rather than 11 before. Of course, we 
observe an excellent energy and angular momentum preservations: the error for 
both quantities is lower than the prescribed tolerance for all parareal iterations 
k > 1. 

In summary, this algorithm, that projects the trajectory on the manifold A4 2 
where both energy and angular momentum are constant, yields results that are 
slightly more accurate than those obtained with Algorithm [2j where the trajec- 
tory is projected on the manifold M. where only the energy is constant. How- 
ever, the behavior of the trajectory is not much improved. From this prototyp- 
ical analysis, it does not seem worth to project the trajectory on the manifold 
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Figure 12: Errors on the energy for Kepler problem, obtained by the parareal 
method ( 25 ]- ( 28 1 with projection on the manifold M.2 of constant energy and 
angular momentum (St = 1(T 4 , dT = 0.01, AT = 0.2). 
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kepler: dt=1 .Oe-4, dT=1 .Oe-2. DT=2.0e-1 , PS + ProjHL 





kepler dt=1 .Oe-4, dT=1 .Oe-2, DT=2.0e-1 , PS + ProjHL 




Figure 13: Errors on the angular momentum for Kepler problem, obtained by 

0.01, AT = 0.2). 



the parareal method (25)-(28) with projection on the manifold M.2 of constant 

10" 4 , dT ■■ 



energy and angular momentum [5t 
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where several invariants (besides the energy) are constant, especially if we bear 
in mind that the determination of these other invariants might be a difficult 
task in the general case. 



6 Symmetric parareal algorithm with symmet- 
ric projection 

We have introduced in Section [5] a parareal algorithm using some projection 
step. Of course, just as the original parareal algorithm Q, this algorithm is 



neither symplectic nor symmetric (see Section 3.2) 



Yet, it is well known that it is possible to couple a symmetric algorithm with an 
appropriate projection step, while still keeping the algorithm symmetric. We 
first briefly recall this idea before using it in our context. 



6.1 Symmetric projection 

Consider the constant energy manifold 

M = {y=(q, P )eR 2d ; H(q,p) = H } , 
and assume that we have at hand a symmetric algorithm with numerical flow 



$at- Consider now the following algorithm (see Fig. 15 for a schematic repre- 
sentation). Assume that !/„ € jM. The approximation y n +i is next defined by 
the set of equations 

= y n + fJ,S7H(y n ), 

= *AT(y„), (29) 
= Vn+l + y^H{y n+ {), 

with (jl chosen such that y n +i <= M.- Consequently, /i and y n +i satisfy the 
equations 




Vn+l - ^AT(y n + pNH{y n )) - fiVH(y n+1 ) = 0, 

H(y n+ i) = H . 



(30) 



We denote by 7r^ m this algorithm: 



M 

Vn+l = 7r %| m (z/n) if and only if (29 1 is satisfied. 



As the same Lagrange multiplier fi is used in the first and third lines of (|29|), 

syir 



and as f at is symmetric, one can easily see that ^TY 1 is a symmetric algorithm. 



To solve the nonlinear equations (30 1, we first recast them in the form S(y n +\, /i) 
0, with 



SiiUn+i,^) \ _ ( Vn+i - ^AT(y n + nVH(y n )) - fiVH(y n+1 ) 
SziVn+i, AO J V H (^A T (y n + ^H(y n j) + fj,VH(y n+1 )) - H 
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Figure 15: Symmetric projection on the manifold Ai. 



We use a Newton-like algorithm, where the jacobian matrix of S is approximated 
by 

JacSto iiA-f 1 -VH(y n )-VH(y n+1 ) \ 

with y n +i — 1 i , AT{yn+t^H(y n ))+fi\7H(y n+ i). To stop the Newton algorithm, 
we again use the three criteria CI, C2, and C3 presented above, where err is 
here the relative residu, 

err = \Si{y n +i, n)\ \S 2 (y n +i,v)\ 

\Vn+l\ \H \ 



An alternative to the algorithm (29) is to use VH(y n+ i) rather than V 'H(y n+ i) 
in the third line of (29). This yields the following algorithm: 

= y n + n'VH(y n ), 

= *AT(y n ), (32) 
= Vn+i + fiVH(y n+1 ), 




with again /i chosen such that y n +i G -M. The advantage in comparison to (29 1 
is that the nonlinear system to be solved is easier. Indeed, once fj, is identified, 
the approximation y n +i can be obtained in an explicit fashion. The Lagrange 
multiplier /j, now solves the equation S(fi) = 0, with 

S(p) = H [*at (y n + MVi?(y»)) + fiVH (y n + f iVH(y n ))} - H , 

which is a scalar nonlinear equation, to be solved for a scalar unknown. We 
denote by sym this algorithm, that we call the quasi- symmetric projection 
algorithm: 

Vn+i = 7r M Sym (y™) ^ anc ^ on ^ ^ 1 32) is satisfied. 

Note that this algorithm is not symmetric. It will yet turn out to be inexpensive, 
and to yield interesting results. 

In practice, the nonlinear problem S(jj,) — is solved using a Newton algorithm, 
where the derivative of S is approximated by 

S'(n) « VH{y n+1 ) T {VH{y n ) + VH(y n )). 
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The algorithm is stopped as soon as one of the three criteria CI, C2, C3 is 
satisfied (in CI, err is the relative error on the energy). 



6.2 Application to the parareal setting 

We now wish to use the above two projection procedures, namely the sym- 
metric projection algorithm tt^ 111 and the quasi-symmetric projection algorithm 
7rT^ sym , in the parareal context. To this end, we return to the formalism used 
in Section P~T1 



We first choose a number K of iterations for the parareal algorithm, and set 

U ■= (u° u 1 ■■■ u K ) 



The symmetric parareal algorithm is given by (171, that we here write as 

U n +i = V%?(U n ), 

where the map ^P^™ is symmetric in the classical sense. We next introduce K 
energies, namely, for any 1 < k < K, 

H k (U n ) :=H(u k n ) 

and we consider a symmetric algorithm, based on the symmetric map 

and on a symmetric projection on the manifold where all the energies H k are 

preserved. Following (29), this algorithm writes 

K 

U n = U n + ^2(i k VH k {U n ) 

k=l 

Un+1 = Virion) 
K 



U n +1 = Un+1 + ^ VkVHk(U n +l), 



k=l 

where the Lagrange multipliers \i k are such that H k (U n +i) = Ho f° r anv 1 < 
k < K. 

Since \7H k (U n ) = (0, . . . , 0, VH {u n ), 0, . . . , 0), the above equations read 
u n = u° n and u n =u k n + {i k VH(u n ) for any 1 < k < K, 

U n+1 = ^(Un), 

u n +i = u n+1 and u n+1 = u n+1 + ^ k VH{u k n+1 ) for any 1 < k < K. 

This yields the following algorithm, that we call symmetric parareal algorithm 
with symmetric projection in the sequel: 

Algorithm 3 Let uq be the initial condition. 
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1. Initialization: set vS. = uq, and sequentially compute \ u° /0 f 

I n+1 /^J 0<n<N-l 

and {<+l} <n<N-l aS 

K+l/2 = Q-AT/2 < y U n), u n+l = GaT/2 («° + l/ 2 ) • 



-Set ui , i /0 = u 



-AT/2 

1/2 — "n+1/2- 



U. Assume that, for some k > 0. i/ie sequence < u k , , , > is known. 

I n ~*~'-i z ) 0<n<N-l 

Compute the sequences {u k+1 \„^ ^„ and < u k+ , \ ,„ > a£ iteration 

>0<n<N { "+ 1 / 2 Jo<n<Ar-l 

fe + 1 &j/ £/ie following steps: 

(a) For allO <n< N—l, compute in parallel T -at/2^\ 1 +i/2) > -^at/2(^ + 
G-AT/2(u* +1/2 ) and 0at/2(w* +1/2 ); 

(c) For < n < N — 1, sequentially compute , 2 and u„Jx a,s 
■ fe+1 - - fc+1 +/i fc+1 Vi J « +1 ) 



w„ = u;, 

- -7 r -AT/2(^ + i/ 2 ) + 5-AT/a(«n+l/a) 
J£l = e?AT/ 2 (^+i /2 ) + ^AT/a^+l/a) - ^AT/2(^ +1/2 ) 



u n+l/2 — y -AT/2 



<\\ = u k + 1 1+f i k+1 WH(u k +\) 
where fik+i is such that = Hq- 

We observe that all the expensive computations (involving the fine propagator) 
can again be performed in parallel. Like in the symmetric parareal algorithm, 
the map that defines u!^\\ from it^ +1 is not symmetric in the classical sense, 
but the map that defines (u^ +1 , u 1 n+1 , • • • , u„ +1 ) from (tt° , u\, • • • , u^) is sym- 
metric. 



Remark 7 Instead of using the symmetric projection scheme (29), we can use 



the quasi-symmetric projection scheme ( 32 I . In the above algorithm, it amounts 
to replacing the last line in 2c), namely 



,k+i 

l n+l 



-rfc+1 



by 



,k+l 
l n+l 



where again fik+i is such that H(u n ~t\) = Hq. This algorithm is called symmet- 
ric parareal algorithm with quasi-symmetric projection. 
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6.3 Evaluation of the complexity of Algorithm [3] 



Under the same hypothesis as in Section |3.3| and the same implementation 
following [3] , we obtain that the complexity of Algorithm [3] is of the order of 



C, 



SP/sym-proj 



SP/sym— proj 



AT 

~5t' 



where Kgp /sym— proj is the number of parareal iterations. The speed-up is hence 
T 

equal to — -— , which is the number of processors divided by the 

SP/sym — proj 

number of iterations. A similar analysis holds for the symmetric parareal scheme 
with quasi- symmetric projection. 



6.4 Numerical results for the Kepler problem 

We first consider Algorithm|3j namely the symmetric projection algorithm, with 
a projection on the constant energy manifold. The iterative projection procedure 
is stopped using the three criteria described ab ove, w ith the parameters tol = 

and 18 for the errors ( 12 1 



~ 7 and N^ K = 2. Results arc shown on Figs. 



16 



17 



10 

on the trajectory, the relative errors (111 on the energy preservation and the 



relative errors ( 27 ) on the angular momentum preservation, respectively. 



We observe that, after only 5 iterations, the errors on the trajectory are com- 
parable to the errors on the trajectory of the fine scheme used sequentially on 
the whole interval [0, 10 4 ]. As expected, the energy is extremely well preserved 
at any parareal iteration k > 1, although, as above, we limit the number of 
iterations to solve the nonlinear projection step. Otherwise stated, convergence 
in this nonlinear procedure is reached for a very small number of iterations. We 
also observe that, for all parareal iterations k > 7, the angular momentum is 
preserved with a relative accuracy of 5 x 10 -4 over the complete time range. 



These results are better than those reported in Section 5.3 for Algorithm [2] 
(namely the standard parareal algorithm coupled with a standard, not sym- 
metrized, projection step). First, trajectory convergence is obtained after 5 



parareal iterations, rather than 11 (compare Figs. 10 and 16 1. Second, numeri 



cal results illustrate the fact that the nonlinear projection step on the constant 
energy manifold converges more rapidly. With our newly constructed Algo- 
rithm [3j the relative error on the energy preservation is smaller than the toler- 
ance 10~ 7 for any parareal iteration k > 1, whereas it is the case only for k > 6 
for Algorithm [2] (compare Figs. [9] and 17). Similarly, the angular momentum is 
better preserved, at any parareal iteration (compare Figs. 11 and 18). 



We next consider the quasi-symmetric projection algorithm, with projection on 
the same manifold, with the same parameters. Results are shown on Figs 
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20 and 21 for the errors (12) on the trajectory, the relative errors (11) on 
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kepler: di=1 .Oe-4, dT=1 .Og-2, DT=2.0e-1 , SPS + SymmProjH 




1 10 100 1000 10000 



kepler dt=1 .Oe-4, dT=1 .Oe-2, DT=2.0e-1 , SPS + SymmProjH 




1 10 100 1000 10000 



;ure 17: Errors on the energy for Kepler problem obtained by Algorithm 
= lO" 4 , dT = 0.01, AT = 0.2). 
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Figure 18: Errors on the angular momentum for Kepler problem obtained by 
Algorithm |] (St = 10~ 4 , dT = 0.01, AT = 0.2). 
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the energy preservation and the relative errors (271 on the angular momentum 
preservation, respectively. 

We observe comparable performances in terms of trajectory accuracy and an- 
gular momentum preservation as with Algorithm [3] (which uses a symmetric, 
rather than quasi-symmetric, projection). However, the preservation of energy 
is not as good: the error on the energy is smaller than the tolerance requested by 
criteria CI on the complete time range [0, T] only for parareal iterations k > 8, 
whereas it is the case for any k > 1 when the symmetric projection is employed 



(compare Figs. 17 and 20 1 




kepler dt=1 ,0e-4. dT=1 .Oe-2, DT=2.0e-1 , SPS + OuasiSymmProjH 



5 0.0001 - 




Figure 19: Errors on the trajectory for Kepler problem obtained by the symmet- 
ric parareal method with quasi-symmetric projection onto the constant energy 
manifold (St = 1CT 4 , dT = 0.01, AT = 0.2). 
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Figure 20: Errors on the energy for Kepler problem obtained by the symmet- 
ric parareal method with quasi-symmetric projection onto the constant energy 
manifold (St = 1CT 4 , dT = 0.01, AT = 0.2). 
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kepler: dt=1 .09-4, dT=1 .Oe-2, DT=2.0e-1 , SPS + QuasiSymmPrajH 




Figure 21: Errors on the angular momentum for Kepler problem obtained 
by the symmetric parareal method with quasi-symmetric projection onto the 
constant energy manifold (St = 1CT 4 , dT = 0.01, AT = 0.2). 
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7 An example in higher dimension: the outer 
solar system 



We now consider a system in a dimension higher than those considered previ- 
ously. We study here the evolution of the outer solar system, which is com- 
posed of 6 three-dimensional particles, representing the Sun and the planets 
from Jupiter to Pluto. The Hamiltonian reads 

H{q,p) = \p T M- 1 p + V{q), q € R 18 , p G R 18 , (33) 



with 



V{q) = - V (34) 



i<j<i<6 



and where M is the diagonal mass matrix M = diag(wi). The values for 
the initial conditions, the gravitational constant G and the masses rrii of the 
particles are taken from [20j Sec. 1.2.4]. By convention, the first particle is the 
Sun, the mass of which is 1000 times as large as the mass of the heaviest planet 
we consider. 

In the sections above, we pointed out that Algorithm [3] yields good results for 
the Kepler problem, and converges in few iterations to the results obtained 
using the fine propagator in a purely sequential manner. We show here that 
(i) this algorithm also behaves well to simulate the outer solar system, whose 
dimensionality is higher, and (ii) that the coarse solver can be driven by a 
simplified dynamics without damaging the good properties of the algorithm. 
The advantage of choosing a simplified coarse integrator is that the complexity 
of the sequential computations using the coarse propagator is now negligible with 
respect to the complexity of the parallel computations using the fine integrator 
on the interval [nAT, (n + 1)AT], in opposition to what was the case for the 
Kepler system. 



The dynamics governed by the Hamiltonian ( 33 1-( 34 ) can indeed be considered 
as a perturbation of the dynamics driven by 



H(q,p) = ~p T M- 1 p+V siuip {q), qeR 18 , p £ M 18 , (35) 



with 



^ II <?!-<£/ II 



EG mi rrij 

2< 3 <6 

In Kiimp, we only take into account the gravitational interaction between the 
Sun (at position qi) and the other planets (at position <jl-, 2 < j < 6), and 



we ignore the interaction between pairs of planets. The fact that (33)-(34) is a 



small perturbation to (35)- (36 1 owes to the huge discrepancy between the mass 
mi of Sun and the masses of the planets. The latter system is less expensive to 
simulate, since V^ mp is a sum of 5 terms, whereas V is a sum of 15 terms. 
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We set St = 10~ 2 , dT = 50 and AT = 200, and study the algorithms over the 
time range [0, T] with T = 2 x 10 5 . 

Results obtained with Algorithm [3] are shown on Figs. [22] [23] and [24] As above, 
we project on the manifold Ai of constant energy, and we use the three cri- 
teria CI, C2, C3 to stop the nonlinear projection step, with the parameters 
^iter* = ^ and tol = 10~ n (which again corresponds to the error on the energy 
preservation of the fine scheme). 

For k > 8, the energy is preserved up to the tolerance prescribed by criteria 
CI on the whole time range. The first convergence criteria that is fulfilled is 
thus CI (hence convergence in the projection step is achieved with less than 
^iter* = 2 iterations). Note that a larger value of N^ x would be required to 
reach convergence for k < 7. We also obtain a good preservation of the angular 
momentum, when k > 5 (the error is then always smaller than 1 %). At k = 15, 
the trajectory error is comparable to the error made by the fine propagator. 
Note that the error on the trajectory is already small when k > 9 (the absolute 



error (12) is smaller than 0.01, and the relative error is even smaller, as \\q\\ is 



of order 1 or larger, see Fig. 25 1 



Remark 8 We have checked that, if we set = 5, we obtain a better en- 

ergy preservation (the error on the energy is then smaller than the tolerance 
prescribed by criteria CI over the complete time range, at any parareal itera- 
tion k ), and comparable results for the angular momentum and the trajectory. 
However, the speed-up is then a little smaller. 



Algorithm [3] also provides very good qualitative results in terms of trajectory. 
On Fig. [25| we plot the actual trajectories of the different planets obtained 
at parareal iterations k = 1, 5, 10 and 15. Except for k = 1, trajectories 
are qualitatively similar to the ones obtained with a symplectic or symmetric 
integration of the solar system. We also note that, even though the trajectory 
is quantitatively wrong for k = 5, it is qualitatively correct. This is a standard 
observation in geometric integration. In turn, for k > 10, the trajectory is 



quantitatively correct, as we observed from Fig. 24 



We conclude this section by discussing the complexity of the algorithm. Note 
that the coarse solver integrates a simpler system than the fine solver. Its 
complexity, for an equal time step, is hence smaller. We are in the case of 



assumption (14). In addition, the average number of iterations in the nonlinear 

m>Gvin — nroi 1.12. 



projection problem is 



With these choices, we need 

-Ksp/sym-proj = 15 parareal iterations 
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outer soler: dt=1 .Oe-2, dT=5.0e1 , DT=2.0e2, SPS + SymmProjH 




outer soler: dt=1 .Oe-2, dT=5.0e1 , DT=2.0e2„ SPS + SymmProjH 



k=11 
k=13 
k=15 



\ a n 



v,"'.':v: : 



10000 100000 
time 



Figure 22: Errors (11 1 on the energy for the outer solar system obtained using 
Algorithm |] (St = 10" 2 , dT = 50, AT = 200). 
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outer soler: dt=1 .Oe-2. dT=5.0e1 . DT=2.0e2, SPS + SymmPrajH 




outer soler: dt=1 .Oe-2, dT=5.0e1 , DT=2.0e2,, SPS + SymmPrajH 



k=7 

k=8 

k=9 

k=11 
k=13 
k=15 

fine 














; vp» r ' 1 ! 1 : i 
"y 1 1 1 



10000 



Figure 23: Errors (27 1 on the angular momentum first component for the outer 
solar system obtained using Algorithm |] (St = 1CT 2 , dT = 50, AT = 200). 
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Figure 24: Errors ( [12] ) on the trajectory for the outer solar system obtained 
using Algorithm ^(St = 1(T 2 , dT = 50, AT = 200). 
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Figure 25: Trajectory of the outer solar system obtained using Algorithm [3] 
(St = 10~ 2 , dT = 50, AT = 200), at iteration k = 1 (top left), k = 5 (top 
right), fc = 10 (bottom left), k = 15 (bottom right). 
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to reach convergence, namely to obtain an error on the trajectory as small as 
that obtained with the fine solver used sequentially over the complete time range 
(see Fig. [24}. 

Using our algorithm, and assuming that we have 1000 processors at our disposal, 
we thus obtain a speed-up of the order of 1000/15 ~ 66: for a similar accuracy, 
the computational time is 66 times smaller. 



8 Conclusion 

The parareal algorithm allows to take benefit from a parallel architecture to 
speed up the numerical integration of systems of ODEs or time-dependent PDEs. 
It consists in the iterative use of a coarse propagator and a fine one. 

In case the system is Hamiltonian, discretization schemes that preserve some 
geometrical features of the underlying Hamiltonian dynamics are known to yield 
the best results, in terms of preservation of energy, but also preservation of other 
first invariants, trajectory accuracy, . . . These properties are not compatible with 
the plain parareal algorithm, which has some issues when used for extremely 
long simulations of Hamiltonian systems. 

We have identified in this article two ingredients that, when combined, yield 
an algorithm that is better adapted to the Hamiltonian context. The first 
ingredient consists in making the scheme symmetric. The idea is simple, but the 
adaptation to the parareal context is not so simple, because the original parareal 
algorithm does not naturally write as a one-step method. We have shown here 
how to obtain a symmetric algorithm that retains the parallel features of the 
original, plain, parareal algorithm. This is yet not sufficient to provide an 
improvement on the geometric properties of the resulting scheme and thus on 
the long time simulations. The second ingredient consists in projecting the 
solution on the constant energy manifold. Here, this projection is not realized 
at every time step but only at the coarse level. This may be the reason why 
such a projection, which is generally not recommended as a viable solution for 
a general solver, appears here to be the right complement (at least in the form 
of the symmetric or quasi symmetric projector) to the symmetrization of the 
parareal algorithm. Note also that the projection is inexpensive, as only one or 
two iterations are needed to solve this nonlinear procedure. 

The conjunction of both ingredients above, in a case where the fine and coarse 
solvers have good geometric properties, is shown here to be the best choice. The 
symmetrized parareal algorithm has, as expected, good geometric features, apart 
from the resonance problem. The projection at the end of each propagation 
interval [T n , T n+ {\ takes care of this resonance problem. After a limited number 
of parareal iterations, we obtain a trajectory which is as accurate as the one 
obtained using the sequential fine solver over the whole time range [0,T]. In 
turn, this convergence stabilizes the values of the invariants, besides the energy. 
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This allows to achieve a substantial speed-up. Let us recall that the full efficiency 
offered by the parallel architecture with the parareal in time algorithm can very 
scarcely be obtained on systems of ODE's. In order to get a full speed-up, the 
parareal algorithm (in its symmetric version or not) has to be combined with 
other iterative procedures (such as domain decomposition approaches for PDEs, 
see e.g. [3"2"]h 

Finally, we emphasize that, in this article, in addition to energy conservation, 
we have chosen to illustrate the quality of the integrators by the accuracy with 
which the trajectories are computed. In many examples of Hamiltonian systems, 
for even longer simulations, the trajectories cannot be well approximated and 
the interest is more in e.g. averages of some quantities along the trajectories. 
The problem should then be considered from a different viewpoint as the one 
adopted here. 

Acknowledgements 

The work of Xiaoying Dai and Yvon Maday is supported in part by the Agence 
Nationale de la Recherche, under the grant ANR-06-CIS6-007 (PITAC). The 
work of Claude Le Bris and Frederic Legoll is supported by INRIA under the 
grant "Action de Recherche Collaborative" HYBRID and by the Agence Na- 
tionale de la Recherche, under the grant ANR-09-BLAN-0216-01 (MEGAS). 

References 

[1] H.C. Andersen, Rattle: a "velocity" version of the Shake algorithm for 
molecular dynamics calculations, J. Comput. Phys. 52, 24-34, 1983. 

[2] L. Baffico, S. Bernard, Y. Maday, G. Turinici, and G. Zerah, 
Parallel in time molecular dynamics simulations, Phys. Rev. E 66, 057701, 
2002. 

[3] G. Bal, On the convergence and the stability of the parareal algorithm to 
solve partial differential equations, in Domain decomposition methods in 
science and engineering, Lect. Notes Comput. Sci. Eng. 40, R. Kornhubcr, 
R. Hoppe, J. Periaux, O. Pironneau, O. Widlund and J. Xu eds., 425-432, 
Springer Verlag, 2005. 

[4] G. Bal and Y. Maday, A parareal time discretization for nonlinear 
PDE's with application to the pricing of an American put, in Recent de- 
velopments in domain decomposition methods, Lect. Notes Comput. Sci. 
Eng. 23, L.F. Pavarino and A. Toselli eds., 189-202, Springer Verlag, 2002. 

[5] G. Bal and Q. Wu, Symplectic parareal, in Domain decomposition meth- 
ods in science and engineering, Lect. Notes Comput. Sci. Eng. 60, U. 



53 



Langer, M. Discacciati, D.E. Keyes, O.B. Widlund and W. Zulehncr cds., 
401-408, Springer Verlag, 2008. 

[6] D.B. Batchelor, personnal communication. 

[7] A. Bellen AND M. Zennaro, Parallel algorithms for initial value prob- 
lems for nonlinear vector difference and differential equations, J. Comput. 
Appl. Math. 25, 341-350, 1989. 

[8] G. Benettin and A. Giorgilli, On the Hamiltonian interpolation of 
near to the identity symplectic mappings with application to symplectic in- 
tegration algorithms, J. Stat. Phys. 74, 1117-1143, 1994. 

[9] K. BuRRAGE, Parallel and sequential methods for ordinary differential 
equations, Numerical Mathematics and Scientific Computation, Oxford Sci- 
ence Publications, The Clarendon Press, Oxford University Press, New 
York, 1995. 

[10] K. BuRRAGE, Parallel methods for ODEs, Advances in Computational 
Mathematics 7(1-2), 1-3, 1997. 

[11] P. Chartier and B. Philippe, A parallel shooting technique for solving 
dissipative ODE's, Computing 51(3-4), 209-236, 1993. 

[12] C. Farhat, J. Cortial, C. Dastillung, and H. Bavestrello, Time- 
parallel implicit integrators for the near-real-time prediction of linear struc- 
tural dynamic responses, Int. J. Numer. Meth. Engng. 67(5), 697-724, 2006. 

[13] P. Fischer, F. Hecht, and Y. Maday, A parareal in time semi implicit 
approximation of the Navier Stokes equations, in Domain decomposition 
methods in science and engineering, Lect. Notes Comput. Sci. Eng. 40, R. 
Kornhuber, R. Hoppe, J. Periaux, O. Pironncau, O. Widlund and J. Xu 
eds., 433-440, Springer Verlag, 2005. 

[14] D. Frenkel and B. Smit, Understanding molecular simulation, from 
algorithms to applications, 2nd cd., Academic Press, 2002. 

[15] M. Gander and S. Vandewalle, On the superlinear and linear con- 
vergence of the parareal algorithm. Proceedings of the 16th International 
Conference on Domain Decomposition Methods, 2005. 

[16] M. Gander and S. Vandewalle, Analysis of the parareal time-parallel 
time-integration method, SIAM J. Sci. Comput. 29(2), 556-578, 2007. 

[17] I. Garrido, B. Lee, G.E. Fladmark, and M.S. Espedal, Conver- 
gent iterative schemes for time parallelization, Mathematics of Computa- 
tion 75(255), 1403-1428, 2006. 

[18] W. HACKBUSCH, Parabolic multigrid methods, Computing methods in 
applied sciences and engineering VI (Versailles, 1983), 189-197, North- 
Holland, Amsterdam, 1984. 



54 



E. Hairer and C. Lubich, The life span of backward error analysis for 
numerical integrators, Numer. Math. 76, 441-462, 1997. 

E. Hairer, C. Lubich, and G. Wanner, Geometric numerical inte- 
gration: structure-preserving algorithms for ordinary differential equations, 
Springer Series in Computational Mathematics 31, 2002. 

P. Joly, Numerical methods for elastic wave propagation, in Waves in 
nonlinear pre-stressed materials, M. Destrade and G. Saccomandi eds., 181— 
281, Springer- Verlag, 2007. 

P. Joly, The mathematical model for elastic wave propagation. Effective 
computational methods for wave propagation, in Numer. Insights, 247-266, 
Chapman & Hall/CRC, 2008. 

J. Laskar, A numerical experiment on the chaotic behavior of the Solar 
system, Nature 338, 237-238, 1989. 

J. Laskar, Chaotic diffusion in the Solar system, Icarus 196, 1-15, 2008. 

B. Leimkuhler and S. Reich, Simulating Hamiltonian dynamics, Cam- 
bridge University Press, 2004. 

B.J. Leimkuhler and R.D. Skeel, Symplectic numerical integrators in 
constrained Hamiltonian systems, J. Comput. Phys. 112, 117-125, 1994. 

E. Lelarasmee, A.E. Ruehli, and A.L. Sangiovanni-Vincentelli, 
The waveform relaxation method for time-domain analysis of large scale 
integrated circuits, IEEE Trans, on CAD of IC and Syst. 1, 131-145, 1982. 

J.-L. Lions, Y. Maday, and G. Turinici, A parareal in time discretiza- 
tion ofPDE's, C. R. Acad. Sci. Paris, Serie I 332, 661-668, 2001. 

Y. Maday, The parareal in time algorithm, in Substructuring Tech- 
niques and Domain Decomposition Methods, F. Magoulcs ed., Chap- 
ter 2, pp. 19-44, Saxe-Coburg Publications, Stirlingshire, UK, 2010. 
doi:10.4203/cscts.24.2 

Y. Maday and G. Turinici, A parareal in time procedure for the control 
of partial differential equations, C. R. Acad. Sci. Paris, Serie I 335, 387-392, 
2002. 

Y. Maday and G. Turinici, A parallel in time approach for quantum 
control: the parareal algorithm, Int. J. Quant. Chem. 93, 223-228, 2003. 

Y. Maday and G. Turinici, The parareal in time iterative solver: a fur- 
ther direction to parallel implementation, in Domain decomposition meth- 
ods in science and engineering, Lcct. Notes Comput. Sci. Eng. 40, R. Ko- 
rnhuber, R. Hoppe, J. Periaux, O. Pironneau, O. Widlund and J. Xu eds., 
441-448, Springer Verlag, 2005. 



55 



[33] A. Quarteroni and A. Valli, Domain decomposition methods for par- 
tial differential equations, Numerical Mathematics and Scientific Computa- 
tion, Oxford Science Publications, The Clarendon Press, Oxford University 
Press, New York, 1999. 

[34] S. Reich, Backward error analysis for numerical integrators, SIAM J. Nu- 
mer. Anal. 36, 1549-1570, 1999. 

[35] J. -P. Ryckaert, G. Ciccotti, and H.J.C. Berendsen, Numerical in- 
tegration of the cartesian equations of motion of a system with constraints: 
molecular dynamics of n-alkanes, J. Comput. Phys. 23, 327-341, 1977. 

[36] P. Saha, J. Stadel and S. Tremaine, A parallel integration method for 
Solar system dynamics, Astron. J. 114(1), 409-414, 1997. 

[37] P. Saha and S. Tremaine, Symplectic integrators for solar system dy- 
namics, Astron. J. 104, 1633-1640, 1992. 

[38] J.M. Sanz-Serna and M.P. Calvo, Numerical Hamiltonian problems, 
Chapman & Hall, 1994. 

[39] G.A. Staff and E.M. R0NQUIST, Stability of the parareal algorithm, in 
Domain decomposition methods in science and engineering, Lect. Notes 
Comput. Sci. Eng. 40, R. Kornhubcr, R. Hoppe, J. Periaux, O. Pironneau, 
O. Widlund and J. Xu eds., 449-456, Springer Verlag, 2005. 

[40] A. Toselli and O. Widlund, Domain decomposition methods- 
algorithms and theory, Springer Series in Computational Mathematics 34, 
Springer Verlag, 2005. 



56 



